function [xh, Sx, oNoise] = ukf_observation(state, Sstate, oNoiseObject, obs, U2, modelObject)

Xdim  = modelObject.statedim;                               % extract state dimension
Vdim  = modelObject.Vdim;                                   % extract process noise dimension
Ndim  = oNoiseObject.dim;                                        % extract observation noise dimension

%Xh is the state
xh  = state;
Sx_ = Sstate;

% Get and calculate CDKF scaling parameters and sigma point weights
h = modelObject.scaleFactor;
hh = h^2;

W1 = [(hh - Xdim - Vdim)/hh   1/(2*hh);                  % sigma-point weights set 1
    1/(2*h)                sqrt(hh-1)/(2*hh)];

W2      = W1;
W2(1,1) = (hh - Xdim - Ndim)/hh ;                        % sigma-point weights set 2


Zeros_Xdim_X_Ndim = zeros(Xdim,Ndim);
Zeros_Ndim_X_Xdim = zeros(Ndim,Xdim);

nsp2   = 2*(Xdim+Ndim) + 1;          % number of sigma points (second set)

Sn = oNoiseObject.cov;         % matrix square root of measurement noise covariance
%------------------------------------------------------
% MEASUREMENT UPDATE

%build up our sigma point set for the observation matrix
Z  = cvecrep([xh; oNoiseObject.mu] ,nsp2);
Sz = [Sx_ Zeros_Xdim_X_Ndim; Zeros_Ndim_X_Xdim Sn];
hSz = h*Sz;
hSzM = [hSz -hSz];
Z(:,2:nsp2) = Z(:,2:nsp2) + hSzM;

%propagate sigma points through observation model
observationState1 = Z(1:Xdim,:); %bascially the current state
observationStateSP = Z(Xdim+1:Xdim+Ndim,:);
Y_ = oNoiseObject.hfun(modelObject,oNoiseObject, observationState1, observationStateSP, U2);

%compose observation covariance (3.231-3.233)
yh_ = W2(1,1)*Y_(:,1) + W2(1,2)*sum(Y_(:,2:nsp2),2);
C = W2(2,1) * ( Y_(:,2:Xdim+Ndim+1) - Y_(:,Xdim+Ndim+2:nsp2));
D = W2(2,2) * ( Y_(:,2:Xdim+Ndim+1) + Y_(:,Xdim+Ndim+2:nsp2) - cvecrep(2*Y_(:,1),Xdim+Ndim));

[~,Sy] = qr([C D]',0);
Sy = Sy';

Syx1 = C(:,1:Xdim);
Syw1 = C(:,Xdim+1:end);

%Calculate the corss coveriance between measurment and observation
Pxy = Sx_*Syx1';

%kalman gain
KG = (Pxy / Sy') / Sy;

%innovation
inov = obs - yh_;

%Kalman update
xh= xh + KG*inov;

%qr update
[~,Sx] = qr([Sx_-KG*Syx1 KG*Syw1 KG*D]',0);
Sx=Sx';


